package model.doe;

import java.util.ArrayList;
import java.util.Date;
import java.util.List;
import java.util.Random;

public class RomeroDickeyGrating extends AbstractDoe{

	private double maxG, startPhase, endPhase;
	private int n, gStep;
	private List<Double> maxAlfas = new ArrayList<Double>();
	private List<Double> currentAlfa = new ArrayList<Double>();
	private AbstractDoe modan;
	private double nsl;
	private double[] phase;
	private double[][] phase2D;
	private double[] xPhase;
	double neededEff;
	int maxIteration;

	public RomeroDickeyGrating(int n, int gStep, double startPhase, double endPhase, AbstractDoe modan) {
		super();
		this.modan=modan;
		this.n = n;
		this.gStep = gStep;
		this.startPhase = startPhase;
		this.endPhase = endPhase;
		maxG = -150000;
		maxAlfas.add(0, 0d);
		currentAlfa.add(0, 0d);
		for (int i = 1; i <= n; i++) {
			maxAlfas.add(i, 0d);
			currentAlfa.add(i, 0d);
		}
	}

	public void buildGratinfWithOptimalParams(double neededEff, int maxIteration)
	{
		List<Double> currentMaxAlfas = new ArrayList<Double>();
		for (int i = 0; i <= n; i++) {
			currentMaxAlfas.add(i, 0d);
		}
		
		int iter = 0;
		double currentEff = -10;
		Random rand = new Random();
		double maxD = 0;
		double currentMaxG = 0;
		double currentd = 0;
		while(iter<maxIteration)
		{
			maxAlfas.set(0, 0d);
			currentAlfa.set(0, 0d);
			currentd = rand.nextDouble()*100000;
			for (int i = 1; i <= n; i++) {
				maxAlfas.set(i, 0d);
				currentAlfa.set(i, currentd);
			}
			computeAlfa();
			computeEff();
			if(currentEff<getEff())
			{
				maxD = currentd;
				currentEff = getEff();
				for (int i = 0; i <= n; i++) {
					currentMaxAlfas.set(i, maxAlfas.get(i).doubleValue());
				}
				currentMaxG = maxG;
			}
			if(currentEff>=neededEff)
			{
				break;
			}
			System.out.println("currentEff " + currentEff + " getEff "+ getEff() + "  iteration " + iter);
			iter++;
		}
		maxAlfas.set(0, currentMaxAlfas.get(0));
		currentAlfa.set(0, 0d);
		for (int i = 1; i <= n; i++) {
			maxAlfas.set(i, currentMaxAlfas.get(i));
			currentAlfa.set(i, maxD);
		}
		maxG = currentMaxG;
		computeEff();
		computePhaseFunction();
	}
	
	private double optimizationFunction(List<Double> alfas, double x) {
		double plsx = 1;
		double qlsx = 0;

		for (int k = 1; k <= n; k++) {
			plsx += 2 * Math.cos(alfas.get(k)) * Math.cos(k * x);
			qlsx += 2 * Math.sin(alfas.get(k)) * Math.cos(k * x);
		}
		return Math.sqrt(plsx * plsx + qlsx * qlsx);
	}

	private double modelS(double x) {
		double plsx = 1;
		double qlsx = 0;

		for (int k = 1; k <= n; k++) {
			plsx += 2 * Math.cos(maxAlfas.get(k)) * Math.cos(k * x);
			qlsx += 2 * Math.sin(maxAlfas.get(k)) * Math.cos(k * x);
		}
		return plsx * plsx + qlsx * qlsx;
	}

	private void computePhaseFunction() {
		try {
			phase = new double[gStep];
			xPhase = new double[gStep];
			for (int i = 0; i < gStep; i++) {
				double plsx = 1;
				double qlsx = 0;
				xPhase[i] = -Math.PI + ((double) i) * 2 * Math.PI
						/ ((double) gStep);
				for (int k = 1; k <= n; k++) {
					plsx += 2d * Math.cos(maxAlfas.get(k))
							* Math.cos(k * xPhase[i]);
					qlsx += 2d * Math.sin(maxAlfas.get(k))
							* Math.cos(k * xPhase[i]);
				}
				phase[i] = Math.atan2(qlsx, plsx);
			}
			for (int i = 1; i < gStep; i++) {
				phase[i]=phase[i-1] + ((phase[i]-phase[i-1])%(2*Math.PI));
				if((phase[i-1]-phase[i])/(xPhase[i-1]-xPhase[i])>=100)
				{
				phase[i]=phase[i] - 2*Math.PI;
				}
				else if(((phase[i-1]-phase[i])/(xPhase[i-1]-xPhase[i]))<=-100)
				{
					phase[i]=phase[i] + 2*Math.PI;
				}
			}
		} catch (Exception ex) {
			ex.printStackTrace();
		}
	}

	private void computeEff() {
		double s = integralS(maxAlfas);
		nsl = maxG * maxG / s;
	}

	private double integralG(List<Double> alfas) {
		double g = 0;
		double x;
		double gx = optimizationFunction(alfas, -Math.PI);
		double gxx = optimizationFunction(alfas, -Math.PI);
		for (int i = 0; i < gStep; i++) {
			double dx = 2 * Math.PI / gStep;
			x = -Math.PI + i * dx;
			gxx = gx;
			gx = optimizationFunction(alfas, x);
			g += (gx + gxx) / 2 * dx;
		}
		return 1 / (2 * Math.PI) * g;
	}

	private double integralS(List<Double> alfas) {
		double s = 0;
		double x;
		double sx = optimizationFunction(alfas, -Math.PI);
		double sxx = optimizationFunction(alfas, -Math.PI);
		for (int i = 0; i < gStep; i++) {
			double dx = 2 * Math.PI / gStep;
			x = -Math.PI + i * dx;
			sxx = sx;
			sx = modelS(x);
			s += (sx + sxx) / 2 * dx;
		}
		return 1 / (2 * Math.PI) * s;
	}

	public void gratingComputing() {
		computeAlfa();
		computePhaseFunction();
		computeEff();
		System.out.println("eff = "+getEff()*100);
	}

	private void computeAlfa() {
		gradientCoordinatesAlfaIteration();
		//simpleCoordinatesAlfaIteration();
	}

	private void simpleCoordinatesAlfaIteration() {
		System.out.println("START SIMPLE");
		Date d1 = new Date();
		maxG = 0;
		double g1, g2;
		double step = 10;
		boolean needFind = true;
		
		int fail;
		while (needFind) {
			fail = 0;
			List<Double> alfas = new ArrayList<Double>();
			for (int i = 0; i <= n; i++) {
				alfas.add(i, currentAlfa.get(i));
			}
			g1 = integralG(alfas);
			double g11;
			for (int i = 1; i <= n; i++) {
				double old = alfas.get(i);
				g11 = integralG(alfas);
				alfas.set(i, alfas.get(i) + step);
				if (integralG(alfas) < g11) {
					alfas.set(i, old-step);
					if (integralG(alfas) < g11) {
					alfas.set(i,old);
					fail++;
					}
				}
			}
			for (int i = 0; i < alfas.size(); i++) {
				currentAlfa.set(i, alfas.get(i));
			}
			double g3 = integralG(alfas);
			if ((Math.abs(g1 - g3) < 0.00000001 && fail != n) || step < 0.00000000000000001) {
				maxG = g3;
				saveAlfa(alfas);
				needFind = false;
			} else if (fail == n) {
				step = step * 0.5;
			}
		}
		Date d2 = new Date();
		System.out.println("TOTAL TIME " + (d2.getTime()-d1.getTime()));
	}

	private void gradientCoordinatesAlfaIteration() {
		System.out.println("START GRADIENT");
		Date d1 = new Date();
		maxG = 0;
		double gBeforeIteration, gBeforeFind;
		double step = 100;
		boolean needFind = true;
		double dx = 0.00001;
		int fail;
		List<Double> alfas = new ArrayList<Double>();
		for (int i = 0; i <= n; i++) {
			alfas.add(i, currentAlfa.get(i));
		}
		while (needFind) {
			
			gBeforeFind = integralG(alfas);
			fail = 0;
			for (int i = 1; i <= n; i++) {
				gBeforeIteration = integralG(alfas);
				//готовим alfas для вычисления частной производной
				alfas.set(i, alfas.get(i) + dx);
				double diffG = integralG(alfas);
				alfas.set(i, alfas.get(i) - dx);
				//частная производная по текущей координате
				double dgdx = (diffG-gBeforeIteration)/dx;
				//начнем искать коэффициент
				double old = alfas.get(i);
				alfas.set(i, old + step*dgdx);
				double gAfterIteration = integralG(alfas);
				if (gAfterIteration < gBeforeIteration) {
					fail++;
					alfas.set(i, old);
				}
				//System.out.println("iter  fail " + fail + "  step  " + step +"  dif  "+ (gAfterIteration - gBeforeIteration) );
			}
			
			
			double gAfterFind = integralG(alfas);
			if ((Math.abs(gBeforeFind - gAfterFind) < 0.00000001 && fail!=n) || step <0.00000000000000000001) {
				maxG = gAfterFind;
				saveAlfa(alfas);
				needFind = false;
			} else if(fail==n)
			{
				step=step*0.1;
			}
		}
		for (int i = 0; i < alfas.size(); i++) {
			currentAlfa.set(i, alfas.get(i));
		}
		Date d2 = new Date();
		System.out.println("TOTAL TIME " + (d2.getTime()-d1.getTime()));
	}
	
	private void gradientCoordinatesAlfaIteration2() {
		System.out.println("START GRADIENT");
		Date d1 = new Date();
		maxG = 0;
		double gBeforeIteration, gBeforeFind;
		double step = 10;
		boolean needFind = true;
		double dx = 0.00001;
		int fail;
		List<Double> alfas = new ArrayList<Double>();
		for (int i = 0; i <= n; i++) {
			alfas.add(i, currentAlfa.get(i));
		}
		while (needFind) {
			
			gBeforeFind = integralG(alfas);
			fail = 0;
			for (int i = 1; i <= n; i++) {
				gBeforeIteration = integralG(alfas);
				//готовим alfas для вычисления частной производной
				alfas.set(i, alfas.get(i) + dx);
				double diffG = integralG(alfas);
				alfas.set(i, alfas.get(i) - dx);
				//частная производная по текущей координате
				double dgdx = (diffG-gBeforeIteration)/dx;
				//начнем искать коэффициент
				double old = alfas.get(i);
				double currentStep = step;
				alfas.set(i, old + currentStep*dgdx);
				double gAfterIteration = integralG(alfas);
				while(gAfterIteration < gBeforeIteration && currentStep>0.01)
				{
					currentStep=0.1*currentStep;
					alfas.set(i, old + currentStep*dgdx);
					gAfterIteration = integralG(alfas);
				}
				if (gAfterIteration < gBeforeIteration) {
					fail++;
					alfas.set(i, old);
				}
				//System.out.println("iter  fail " + fail + "  step  " + step +"  dif  "+ (gAfterIteration - gBeforeIteration) );
			}
			
			
			double gAfterFind = integralG(alfas);
			if ((Math.abs(gBeforeFind - gAfterFind) < 0.00000001 && fail!=n) || step <0.00000000000000000001) {
				maxG = gAfterFind;
				saveAlfa(alfas);
				needFind = false;
			} else if(fail==n)
			{
				step=step*0.1;
			}
		}
		for (int i = 0; i < alfas.size(); i++) {
			currentAlfa.set(i, alfas.get(i));
		}
		Date d2 = new Date();
		System.out.println("TOTAL TIME " + (d2.getTime()-d1.getTime()));
	}
	
	private void saveAlfa(List<Double> alfas) {
		for (int i = 0; i < maxAlfas.size(); i++) {
			maxAlfas.set(i, alfas.get(i));
		}
	}

	@Override
	public double run(int i, int j) {
		return 0;
	}
	@Override
	public double [][] getAbsoluteAmplitudeMatrix()
	{
		double[][] matrix = null;
		if(modan!=null && modan.getWidthPointsNumber()==gStep &&modan.getHeightPointsNumber()==gStep)
		{
			matrix = modan.getAbsoluteAmplitudeMatrix();
		}
		else
		{
			matrix = new double[gStep][gStep];
			for (int j = 0; j < gStep; j++) {
				for (int i = 0; i < gStep; i++) {
					matrix[i][j]=1;
				}
			}
		}
		return matrix;
	}
	
	@Override
	public double [][] getPhaseTo2PIMatrix()
	{
		System.out.println("getPhaseTo2PIMatrix");
		double[][] modanPhase = modan.getSignPhaseMatrix();
		phase2D = new double[gStep][gStep];
		double[] currentXPhase = new double[gStep];
		double[][] res = new double[gStep][gStep];
		//собираем фазу делителя
		for (int j = 0; j < gStep; j++) {
			for (int i = 0; i < gStep; i++) {
				double plsx = 1;
				double qlsx = 0;
				currentXPhase[i] = startPhase + ((double) i) * (endPhase - startPhase)
						/ ((double) gStep);
				for (int k = 1; k <= n; k++) {
					plsx += 2d * Math.cos(maxAlfas.get(k)) * Math.cos(k * currentXPhase[i]);
					qlsx += 2d * Math.sin(maxAlfas.get(k)) * Math.cos(k * currentXPhase[i]);
				}
				phase2D[j][i] = Math.atan2(qlsx, plsx);
			}
		}
		//добавляем фазу модана для Фурье
		for (int j = 0; j < gStep; j++) {
			for (int i = 0; i < gStep; i++) {
				res[j][i] = phase2D[j][i] + modanPhase[j][i];
			}
		}
		return res;
	}
	
	public double getEff() {
		return nsl;
	}

	public double[] getPhase() {
		return phase;
	}

	public void setPhase(double[] phase) {
		this.phase = phase;
	}

	public double[] getxPhase() {
		return xPhase;
	}

	public void setxPhase(double[] xPhase) {
		this.xPhase = xPhase;
	}

	public double[][] getPhase2D() {
		return phase2D;
	}

	public void setPhase2D(double[][] phase2d) {
		phase2D = phase2d;
	}

}
